# Figure 9 code

library(RColorBrewer)
library(cshapes)
#set the working directory here
setwd("")

# load the IGO membership array:
load("membership-expanded3.rda")

# load the communities data:
load("clubs of clubs data 2.rda")

# load the IGO statistics:
load("igostats5.rda")

source("ccode2cow.R")
source("report.progress.R")

# first make an array containing the number of shared IGO memberships across dyads 
shared<-array(NA, dim=c(dim(membership)[1], dim(membership)[1], dim(membership)[3]), dimnames=list(dimnames(membership)[[1]], dimnames(membership)[[1]], dimnames(membership)[[3]]))
for (i in 1:dim(shared)[3]){
	shared[,,i]<-membership[,,i] %*% t(membership[,,i])
	diag(shared[,,i])<-NA
	}

# Now crawl through the comm data, calculating the proportion shared for each country-year:
comm$bonding<-NA
for (i in 1:nrow(comm)){
	report.progress(i, nrow(comm))
	cow.current<-comm$cow[i]
	year.current<-comm$year[i]
	comm.current<-comm$community[i]
	memberlist<-comm$cow[comm$year==year.current & comm$community==comm.current]
	totalties<-sum(shared[cow.current, , as.character(year.current)], na.rm=TRUE)
	inwardties<-sum(shared[cow.current, memberlist, as.character(year.current)], na.rm=TRUE)
	comm$totalties[i]<-totalties
	comm$inwardties[i]<-inwardties
	comm$outwardties[i]<-totalties-inwardties
	comm$bonding[i]<-inwardties/totalties
		}

# Do a test check for USA in 1950:
LA1950<-comm$cow[comm$year==1950 & comm$comm.name=="Latin America"]
total<-sum(shared["USA",,"1950"], na.rm=T)	
inward<-sum(shared["USA", LA1950,"1950"], na.rm=T)
outward<-sum(shared["USA", !(dimnames(shared)[[2]] %in% LA1950), "1950"], na.rm=T)
bonding<-inward/total
comm$totalties[comm$cow=="USA" & comm$year==1950]==total # TRUE
comm$inwardties[comm$cow=="USA" & comm$year==1950]==inward # TRUE
comm$bonding[comm$cow=="USA" & comm$year==1950]==bonding # TRUE

# Now build Fig 9:
yearlist<-sort(unique(comm$year))
commlist<-sort(unique(comm$comm.name))
quantiles<-c(0, 0.25, 0.5, 0.75, 1)
b<-array(NA, dim=c(length(yearlist), length(quantiles), length(commlist) ), dimnames=list(year=yearlist, quantiles=quantiles, community=commlist))
for (y in yearlist){
	for (m in commlist) b[as.character(y), , m]<-quantile(comm$bonding[comm$year==y & comm$comm.name==m], quantiles, na.rm=T)
	}
	
pdf("Figure 9.pdf", width=5, height=8)
par(mfrow=c(3,1))
commnames<-c("African", "European/Northern", "Latin American", "South Asian/Other")
for (i in 1:3){
	plot(yearlist, yearlist, type="n", xlab="Year", ylab="Proportion of Intra-Cluster Ties", ylim=range(b, na.rm=T), las=1, bty="n", main=commnames[i])
	#abline(h=seq(0.1, 0.6, 0.1), lwd=0.5, col="grey80")
	for (q in quantiles) lines(yearlist, b[,as.character(q),i], col="grey80")
	lines(yearlist, b[,"0.5",i], lwd=2)
	}
dev.off()